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NOMENCLATURE 

C specific heat (Btu/lb°F) 

2 

h heat transfer coefficient (Btu/hr ft °F) 
k thermal conductivity (Btu/hr ft °F) 

L slab thickness (ft) 

3 

q volumetric heat source (Btu/hr ft ) 

2 2 

q" heat source per unit area (Btu/hr ft or Watts/in ) 
T temperature (°F) 

X space direction (in) 

X dimensionless space direction 

Greek Letters 

2 

a thermal diffusivity (ft /hr) 

At size of time step (sec) 

AT temperature semi-interval across T^ (°F) 

Ax grid spacing in x direction (in) 

AX dimensionless grid spacing in x direction 
A0 dimensionless temperature semi-interval across 6^. 

0 dimensionless temperature 

X latent heat of ice (Btu/ft^) 

p density (Ib/ft^) 

Subscripts 

j layer in the composite body 

1 grid point 

s solid or ice region 



w liquid or water region 
f phase change 

o ambient at lower boundary of composite body 

“ ambient at upper boundary of composite body 

Superscripts 
n time step 

* ice-water "mushy" region 
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I. INTRODUCTION 

The formation of ice on aircraft components, a wing 
section of which is shown in Figure 1, poses tremendous 
difficulties in aircraft operation. Various anti-icing and 
de-icing methods have been investigated and are reported in 
Reference (1). Before these methods are considered, it is 
necessary to differentiate between anti-icing and de-icing. 
The anti-icing principle involves methods to prevent ice 
formation on the blade surface. These methods usually 
require an excessive amount of energy and are not frequently 
applied. The de-icing principle involves shedding the ice 
by heating the surface on which it is formed. In this case 
the energy required is reduced significantly because heat 
is needed only to form a thin water film between the air- 
craft structure and the ice, thereby decreasing the adhe- 
sion strength and allowing the aerodynamic forces to sweep 
away the ice layer. 

At present, for anti-icing and de-icing, thermal energy 
is the most commonly used technique and is obtained in two 
ways : 

(i) Electrical resistance heating elements which are embedded 
just below the surface being heated (as shown in Figure 
1 ) ; 
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(ii) Passing hot, engine compressed bleed air through pas- 
sages below the surface being heated. 

A. TYPES OF DE-ICING AND ANTI-ICING SYSTEMS 

Various de-icing and anti-icing techniques have been 
developed over the years, some of which have been in use for 
a long time and are listed below. 

1. Inflatable Rubber Pneumatic Boots: This pneumati- 

cally operated mechanical de-icing system consists of a boot 
which is made of a flexible rubber-like material, and which 
is slipped over the wing of the aircraft such that the ice 
forms on the boot-surface rather than on the wing structure. 
The boots, when inflated with the cooled engine bleed air, 
break the ice surface, thus allowing the aerodynamic forces 
to blow the ice away. A section of the aircraft wing fitted 
with a boot is shown in Figure 2. This method is relatively 
simple and uses a very small quantity of the cooled engine 
bleed air and does not impose any fuel penalty on the engine. 
However, because the surface of the boot is not as smooth 
as the wing surface, the boot system increases the drag for 
the aircraft. In addition, this pneumatic de-icer system 
requires frequent maintenance and replacement. For the 
pneumatic boots to show an operational weight advantage over 
the thermal de-icing system technique, they must be vir- 
tually drag-free. Modern boots, carefully bonded to the 
wing structure, do minimize the drag but are expensive to 


maintain . 
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2. Chemical Systems: An anti-icing method used to 

prevent the formation of ice on the abrasion shield surface 
of the aircraft involves the use of freezing point depres- 
sants. The chemical depressants are spread in a thin film 
over the abrasion shield, thus lowering the water freezing 
point and preventing the formation of ice. Though relatively 
simple in concept, their usefulness is restricted because 

of the variable external pressure field, which makes it 
difficult to obtain a uniform flow distribution of the de- 
pressants, and which therefore can result in the formation 
of ice on the unprotected areas. Further drawbacks are that 
the systems require frequent resupply and are sensitive to 
clogging of the fluid distribution holes in dusty environ- 
ments. Chemical systems are therefore restricted in their 
use to mainly windshield ice protection. 

3. Thermal De-icing: Thermal de-icing remains the 

most commonly used technique in removing the ice from the 
abrasion shield surface. In this method the ice covered 
regions are cyclically heated in sequence either by electric 
heaters or by hot bleed air from the engine. The thermal 
energy supplied is used to raise thd temperature of the sur- 
face on which the ice is deposited to 32 °F and to melt a 
thin layer of ice. This thin film of water reduces the 
adhesion strength of the ice to the surface and aerodynamic 
forces then sweep the unmelted ice from the surface. Because 
of the cyclical nature of the energy input, the energy 
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requirement is very much less compared to the other techni- 
ques. In addition, the use of this system is not restricted 
due to change in weather conditions and is relatively easy 
to maintain. 

Of the above mentioned systems, the electro- thermal 
de-icing technique is the most commonly used and will be the 
one considered in this study. 

B. DE-ICING PAD CONFIGURATION 

The configuration of the electro-thermal de-icing pad 
is shown in Figure 1. It is essentially a composite body 
consisting of five layers in the case of a point heat source 
and six layers in the presence of a finite heater. The 
heating source is separated from the metal substrate, or 
the aircraft blade, by the inner insulation which usually 
consists of resin impregnated glass cloth. This insulation 
serves to provide electrical insulation for the heating 
element, and also directs the heat towards the ice layer. 

It is advantageous to have a large ratio of inner to outer 
Insulation thickness so that more heat flows toward the 
ice layer, thereby reducing the de-icing time. The heater 
element usually consists either of a woven mat of wires and 
glass fibers or of multiple strips of resistance ribbon. 

In order to protect the de-icer pad from rain erosion 
or sand and stone abrasion, which could be a problem while 
flying at high speeds, an abrasion shield, frequently made 
of stainless steel, is added to the outer insulation. The 
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abrasion shield also serves to diffuse the heat from the 
heater, thus providing more uniform heating and thereby 
reducing cold spots where ice could form above the gaps in 
the heater elements. 

The material of construction of the substrate depends 
on the type of aircraft and is most often an aluminum alloy. 

In this study, only the one-dimensional model of the 
de-icing pad will be investigated. It is assumed that there 
is perfect adhesion between each of the layers and therefore 
no contact resistance to heat flow will be considered in 
this analysis. This study will concern itself with both a 
point heat source and a finite thickness heater, providing 
either constant or time dependent heat output. In each of 
these cases, the effects of the type and thickness of the 
insulation layers and the nature of the blade structure on 
the de-icing time will be observed. In the first section 
of this study, the phase change at the ice-abrasion shield 
interface will not be considered and the ice layer will be 
treated as a single phase. In the latter part, for the 
ice-water phase change, a numerical method which approxi- 
mates the latent heat effect by a large heat capacity over 
a small temperature interval will be applied. 
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II. LITERATURE REVIEW 

From the conducted literature review, it is evident 
that the de-icing problem has either been ignored or that 
information pertaining to it has not been published in the 
open literature. Of the few attempts that have been made 
to solve the specific problem, Wardlaw(2) and Campbell (3) 
applied an analytical approach, while Stallabrass (4) has 
made use of a numerical technique. However, several 
methods have been proposed to solve the transient heat con- 
duction problem in a composite slab, without phase change 
and with different boundary conditions, and these will be 
reviewed below. 

A. ANALYTICAL' TECHNIQUES FOR COMPOSITE BODY HEAT TRANSFER 

The Laplace transformation technique was presented by 
Carslaw and Jaeger (5), but this method becomes increasingly 
tedious to apply to composite bodies with more than two 
layers since it becomes more difficult to obtain the inverse 
Laplace transform. Goodman (6) introduced the method of the 
adjoint solution, which arises from consideration of an 
auxiliary function. However, a disadvantage of the adjoint 
method is that the solution provides just the interfacial 
temperatures and not the temperature profile within each 
slab layer. Further work on the adjoint solution method has 
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been done by Bouchillon ( 7) in which transient cases have 
been considered. The formulated integral equations have 
been reduced to linear equations and then solved by matrix 
inversion techniques. However, as before, only the inter- 
facial temperatures can be calculated. 

The Orthogonal-Expansion technique proposed by Tittle 
(8,9) is another method to solve the boundary-value heat 
conduction problem in multilayer regions. The method is 
basically an extension of the Sturm-Liouville problem to the 
case of a one-dimensional multilayer region. Orthogonal sets 
are constructed from the solution of each of the layers and 
an orthogonality factor, called the discontinuous weighting 
function, is used such that the resulting orthogonal set is 
applicable to the entire composite media. 

Bulavin and Kashcheev (10) used the method of separation 
of variables and of orthogonal expansion of functions over 
a one-dimensional multilayer region to solve the transient 
heat conduction problem involving heat sources in a multi- 
layer region. Campbell (3) applied a similar method in 
solving the de-icer pad problem analytically. 

The disadvantage of using an analytical technique is 
that for each temperature desired, an excessive amount of 
calculations have to be performed. Hence, as the number of 
layers within the body increases, the calculations become 
more tedious. This drawback can be overcome by using numer- 
ical techniques. 
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B. NUMERICAL TECHNIQUES FOR COMPOSITE BODY HEAT TRANSFER 

The most frequently used numerical method for solving 
partial differential equations is the finite-difference 
method. Using this method, the temperature at all the nodal 
points within the composite slab can be calculated at each 
time step. The finite-difference method involves the con- 
struction of a grid within the boundary over which the dif- 
ferential problem is to be solved. At each of the grid 
points the differential operators are replaced by their 
approximate values expressed in terms of functions. This 
substitution reduces the problem to the solution of a set 
of algebraic equations, which is mathematically easier to 
solve. A finite-difference representation of the de-icer 
pad is shown in Figure 3* 

Two of the major considerations in using a finite-dif- 
ference scheme are the establishment of both a convergence 
criteria and a stability criteria. A number of finite- 
differencing methods have been proposed and are discussed 
in depth in References (11) and (12) . In addition. Price 
and Slack (13) have evaluated the accuracy and stability, 
criteria of various finite-differencing methods for the 
heat flow equation with convective boundary conditions. 

The Crank-Nicolson implicit finite-difference scheme 
is unconditionally stable for all time steps, and has been 
used in the present study. For the corresponding explicit 

2 

formulation, there exists a limitation on the ratio At/ (Ax) 
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2 

to be 0 £ At/ (Ax) £ 1/2; no such conditions are present for 
the implicit scheme. However, the choice of the size of the 
time and space steps has a direct effect on the accuracy of 
the solution. The truncation error for the implicit Crank- 
Nicolson scheme is of the order (Ax) for the space differ- 
ential and (At) for the time differential. As stated by 
Dusinberre ( 14) , the accuracy of the solution can be increased 
by initially choosing a very small time step and then sub- 
sequently increasing it. 

Stallabrass (4) employed the explicit finite-difference 
scheme in his de-icer analysis. Though the scope of his 
study covers most of the important aspects for the one- 
dimensional de-icer design, its major drawback is that it 
considers no phase change within the ice layer and limits 
the solution to the time at which the interfacial tempera- 
ture between the ice and the abrasion shield reaches 32°F. 
Also, Stallabrass (4) considers only a point source and a 
finite thickness heater with constant heat output. The 
present study will consider the phase change effect in the 
ice-layer as well as time varying heat sources. 

C, METHODS TO DESCRI.I3E PHASE CHANGE 

Phase change or moving boundary problems have been rela- 
tively difficult to solve because of the non-linear nature 
of the boundary conditions arising from the boundary move- 
ment. Various methods of analysis such as the Integral 
method (15), Success i.ve-Approximation technique (16) and 
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Series Solution (17) have been proposed for analytically 
solving the one-dimensional phase change problem. Mori and 
Araki(18) have reviewed some of the other methods that have 
been proposed. 

Numerical methods to solve the phase change problem 
have been attempted and are described by Meuhlbauer and 
Sunderland ( 19 ) and Rubinstein ( 20) . Most of the numerical 
methods solve the pertinent heat conduction equations and 
determine the temperature distribution in both media, while 
at the same time locating the position of the solid-liquid 
interface by a predictor-corrector technique. However, this 
requires a large number of iterations to locate the solid- 
liquid interface position at any given time. Mastanaiah (21) 
used such an iterative scheme with a two time level implicit 
method for the one-dimensional freezing and melting problem 
with convective boundary conditions and variable thermal 
properties. Lazaridis (22) used another iterative solution 
for the two-dimensional solidification problem with constant 
thermal properties and convective boundary conditions, and 
also constant temperature conditions at the boundaries. 

Crank (23) describes two additional ways to approach the 
moving boundary problem. The first involves the rearranging 
of variables such that the boundary is treated as stationary 
and the problem is transformed into an eigenvalue problem 
with fixed boundaries. However, the equations to be solved 
contain parameters associated with the moving boundary 
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problem for which values have to be determined to satisfy 
the boundary conditions. In the second method, Lagrangian 
interpolation formulae for non-equal intervals are intro- 
duced along with the finite-difference formulae in order to 
follow the movement of the boundary. 

To avoid the problem of locating the interface posi- 
tion as in the above method, a second approach, the method 
of weak solution, often called the Enthalpy Method, has been 
used. In this method, the enthalpy is used as a dependent 
variable along with the temperature. Thus the moving bound- 
ary problem can be solved in a fixed region and no modifi- 
cation is required to satisfy conditions at the moving 
boundary. Much of the numerical work applying the enthalpy 
approach to the phase change (Stefan) problems has been 
done using the finite-difference scheme. Atthey(24) has 
solved the welding problem in one-dimension, which is essen- 
tially a melting problem, using this approach. The conver- 
gence criteria for such a solution has also been clearly 
indicated. The latent heat has, however, been assumed to 
be evolved at the phase change temperature. In practice, 
the latent heat is usually considered to be evolved over a 
small temperature range, AT. 

Goodrich(25) and Bonacina et al. (26) have solved the 
one-dimensional ice-water problem by associating the latent 
heat effect with a finite temperature interval about the 
phase change isotherm. However, it should be noted that 
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within this "mushy" region, the grid spacing has to be sub- 
stantially reduced or else the isotherm may advance in an 
oscillatory fashion and distort the temperature profile. 
Goodrich (25) used the Crank-Nicolson implicit finite-differ- 
ence scheme for formulating the problem and the Gaussian 
Elimination technique for solving the resulting set of equa- 
tions. Bonacina et al. (26) used a three-time level implicit 
scheme for formulating the problem, which was then solved 
as before. The formulation results in three governing equa- 
tions applicable to the three phase regions; solid, "mushy" 
and liquid regions, respectively. The "mushy" region was 
defined over a small finite temperature range, 2AT, about 
the phase change temperature, T^. The phase change initially 
starts occurring at temperature T^ - AT, and the ice becomes 
pure liquid at temperature T^ + AT, where T^ equals 32° F. 

The choice of the temperature interval, 2 AT, depends on the 
physical nature of the problem. For the ice-water system 
considered in Reference (26) , a temperature interval of two 
degrees Kelvin was assumed and good comparison to the ana- 
lytical solution was obtained. In the solid and liquid 
regions the thermal properties were assumed constant. In 
the "mushy" region, the thermal conductivity was assumed to 
vary linearly with temperature and the latent heat effect 
was approximated by a large heat capacity. Solutions to 
phase change problems applying this method agree fairly 
accurately with analytical results and this method will there- 
fore be used in the present study. 
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The enthalpy method has also been employed for solving 
two-dimensional phase change problems. Meyers (27) and 
Shamsunder and Sparrow (28) have described a purely implicit 
two-dimensional finite-difference scheme for solving such 
problems. Shamsunder and Sparrow (28) have also considered 
the effect of various parameters on the solidification rate. 
A finite element approach has been proposed by Comini et al. 
(29) for solution of the Stefan problem in two-dimensions 
with non-linear radiation boundary conditions. 

The complete numerical formulation of the de-icer prob- 
lem is given in the next section. 
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III. NUMERICAL FORMULATION 

A. GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 

In the formulation for the one-dimensional, unsteady 
state, mathematical model for the heat transfer analysis of 
a composite aircraft blade with an ice layer, as shown in 
Figure 1, the following assumptions were made: 

(1) The thermal properties of the material composing each 
layer of the blade are constant; 

(2) Density variations are neglected, as are the effects of 
the volume contraction experienced when the ice melts; 

(3) The individual layers are in perfect contact with each 
other, and there is no additional resistance present at 
the interface; and 

(4) The ambient temperature is constant. 

With the above assumptions the governing equation is 

( 1 ) 

3 t 9 

where j represents the layer in question, p, C, k, T and q 
are the density, heat capacity, thermal conductivity, temp- 
erature and volumetric, heat source of the layer, respec- 

tively, and X and t are the distance coordinate and time 


variable . 
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The total blade is then characterized by: 


j = 1 Blade or Substrate = 0 

j = 2 Lower or Inner Insulation q 2 = 0 

j = 3 Heater q^ = f (t) 

j = 4 Upper or Outer Insulation q^ = 0 

j = 5 Abrasion Shield q^ = 0 


( 2 ) 


For the ice layer (j=6) , the governing equations to be 
applied depend upon the temperature profile within the layer. 
Applying the method discussed at the end of Chapter II, 
References (25,26), the governing equations are: 

For j = 6 , 


Ice ; 


Water: 



,-Jj 


Ice-Water : 


where 


3t 

C‘ = 



Zb.1 2 

A +.^Ti_(VATj) 


In the above equations, and C^, and and and 
k^ are the specific heats, densities and thermal conductivi- 
ties of ice and water, respectively. The latent heat X, is 
assumed to be evolved over the fictitious temperature inter- 
val 2AT (2°F in this study) . The phase change initially 
starts occurring at temperature T^ - AT and the ice becomes 
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pure liquid at temperature T^: + AT, where equals 32° F. 
Within the temperature range - AT and T^ + AT, a single 
phase does not exist and the ice-water system is said to be 
in the "mushy" region. As discussed previously, a large 
heat capacity to account for the latent heat and a linear- 
ized thermal conductivity are used to describe this "mushy" 
region. The thermal properties in the ice and water regions 
are assumed to remain constant. 

The corresponding boundary conditions to be used are: 

(i) The equality of the temperatures and heat fluxes at the 
interfacial points: 




where the subscript I denotes the interface; 

(ii) Convective heat transfer at the lower and upper bound- 
aries : 


k i 


X«0 




(5a) 


where 


■ 

X«L 

= 

Tjl - 

H 

11 


■T-i 


3 

II 

v 

■nl 

Je. 

= r 

• 5-0 r 

T^-^T 


» j = 6 


(5b) 


< X-AT 

■T 




X 
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To formulate the above equations in terms of non-dimen- 
sional temperature and distance, the following definitions 
are made: 

0=X y=A <i'> 

where = the reference temperature (taken to be 32 °F in 

this study) , 

L = the total length of the composite slab, and 

o.^ = the thermal diffusivity of the j " layer. 

Substitution of the above dimensionless quantities into 
equations (1) through (5) yields 


^0 


9 

For j = 6 


4 _ 

t e 9x" 






3 “ lf*»**/5 


( 2 ') 


Ice ; 

30^ 




9 t 

1} 


Water: 

9©i 

o^w» 



dt 

IF 

3 X*- 

Ice-water: 

llcdQi 

- 2. 



3 1 


jkJ 


> e^+t^e 
e.-^e $0: 

^ i ^ 


where 






Js 




I -7^ Q ^ C> 

z 
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At the interfacial points. 


A- 30 L 
3X 


Vfi- 




j 1,>...,5 


X' 


Finally, at the lower and upper boundaries 
at X = 0, j = 1 


* 3X 


X°0 




at X = 1, j = 6 



- 

X=>jL 

K L 
1 

where 




11 

t 

•^ror 


/fej = J 

k* 

:S-ov 


_ . 4,,l (ei|_^-e.) 


Ix«l 


( 4 ') 


(5a’) 


(5b') 


V// • •-♦V/ 

R=l * 

i-ov 6 -Ae ^ a| 6 +a9 


The above equations can now be represented in finite- 
difference form and solved numerically. 

B. NUMERICAL TECHNIQUE 

In the numerical solution of the above partial differ- 
ential equations, a finite-difference scheme was adopted. 
This is accomplished by constructing a system of grid points 
which define a finite number of regularly spaced values of 
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the dependent variable, in this case the temperature, over 
the whole space domain at each time step. As illustrated 
in Figure 3, the X-axis represents the length along the com- 
posite body and the Y-axis represents the time variable. 

The space axis is divided into increments of size AX ^ , where 
the subscript j indicates the layer in question. Within 
each layer, the space increment AX^ is constant, but it may 
vary from one layer to the next. The time axis is divided 
into equal time step intervals. At. However, the time incre- 
ment may be increased as the solution progresses. The index 
i denotes the position of the variable in the grid along the 

space axis, and the superscript n indicates the value of 

n 

the variable at time n. Hence the quantity 0. . represents 

D ^ ^ 

the non-dimensional temperature in the j layer, at posi- 
tion i along the space axis, at time n. 

C. FINITE-DIFFERENCING AND THE METHOD OF SOLUTION 

The finite-difference approximation to a partial deriv- 
ative can be derived using a Taylor Series expansion around 
any grid point. In this study, the Crank-Nicolson Implicit 
scheme was applied in order to maintain stability of the 
solution. The finite-difference approximations for the 
first and second-order space derivatives and first order 
time derivative are second order correct and are given below. 



9X 



+ 0(ax/ 


L 


(6a) 
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sx** 




Mi 

dt 


= 6 j- ;t-»l ^ ^i, t H- Of; L-X+ Oj, L-Vl /l 0 C-1 

+ 0 (AS/ 

-W+l 2 . 

= ~~ Qj-A O(A-t) 

L At 


(6b) 


(6c) 


In the above equations the superscripts represent the 
time levels n and n+1, the subscript j represents the j " 
layer and the subscript i represents the grid point under 
consideration . 

The Crank-Nicolson finite-difference equations are 
obtained by substituting (6a,b,c) into equations (2') 
through (S'). The governing equations in finite-difference 
form reduce to; 






( 7 ) 


For j = 1 5 

where q" is the source per unit area and equals q^'LAX. 

For j = 6, 


x\ Y^+i ^yy^i. 

l-Vl’ 






(8a) 


0j< e^“A9 
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+ ^{^~ -^"J ®i/^ “ 


- 6 ^- 


- _71+i 


+AX:;,- 

-^' eli-i+ (^^--»^'- ©■ , - 

where ” As ^ 

. ■• “• 


v» + )^ . -.'I^^l. 


^ ©;,{.-»+ ©i/tl/^e, -6 ej 

2A8 L 2 *■ J 


C -_>!i 

2Tv<4,A.e 


The value of the temperature at the half level in time 
is obtained from a truncated Taylor Series as follows; 


-Trt + l/j Y» ^ .2 

+ ¥ (Ify + 


The time derivative in the above analog is obtained 
from the equation (3') as 

L'c'heA „ A* ( el^,~ el. ) . a. (el.- eL-) ^ 

6^x)*’ ^ ("Axf / 
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where 





4 - 


r* ^ 

9 ^U.-t-©^,L _ 


;iAe 

L 

1 

A 

4 - 


- >» 1 

[e 7 .i 


<1 

I 

£A© J 


z 

V" A 

A 


t 

■^sCs -h 

/ujCw 


Z A 0 


z 




Substituting (8e) into (8d) yields the finite-difference 
n+Jj 

analog for 0 • , as 

") f J- 


011 =. 6: „ 

ZCHl 



(8f) 


n+ij n+^ 

For 0. , , and 0. . ^ the finite-difference analogs 
1 / i+x ] / i~i 

are obtained in a similar manner and are given below. 


e],- = e” 


^ VITL \ J V (AX/ / 


(8g) 


where 9i]uz+ ^,U, - A ©'^ 

ZA6 £ - - V ^ / 


= 


£ 

•^s ■+* ^ 6j., i-4-i ~t &i,L _ ^0^ _ J 

-/sG+ 


C = A 


/ Tyt^ A© 


-i- 


■Vl +• '/j_ y, 


(8h) 
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Li ' /J 


■0ic-.+ sli-i _ ( flt-z.e'il 

2 A© 

z \ /I 


Z 

2 


Equation (7) is valid for all grid points within each 

^ V*k 

of the j layers except for the ice layer, within which 
equations ( 8a, b, c) are applicable. Equation (8c) was 
obtained using the finite difference formulation suggested 
by Von Rosenberg (12) and the temperatures at the half levels 
in time were calculated applying the method devised by 
Douglas (30) . 

The finite-difference equations for the interfacial 
points, for the two boundary conditions and for the heater 
are discussed below. 

1. Finite-Difference Equations At Interfacial Points 

Let i be the interfacial point between the slab layers 
j and j+1 as shown in Figure 4a. At this point, the temper- 
ature and the heat fluxes are equal and equation (4') is 
valid. Substitution of the finite-difference equation (6a) 
into equation (4') yields; 

• = 0 ^ . 


y\ 



(9a) 
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th 

Since the boundary extends only to the i grid point 
th ^ 

for the j slab, the temperatures 0. . ^ and 6. are 

1 » 1+1 3 , 1+1 

ri+ X 

fictitious. In the same manner, 9.,, . , and 0.,, • , are 

3 + 1 , 1-1 3 + 1 , 1— 1 

fictitious. At the point i, the governing equation (7) 
when applied for the layers j and j+1 yields: 

.. — 2(1+0-^) \0- • -hfi' • - -0‘ . 


_ %uh% ) 4 ', ^ z/i . e]j,, 


(lOa) 




(lOb) 


Eliminating the fictitious temperatures between equa- 
tions (9b, 10a, b) and using equation (9a) yields the following 
applicable finite-difference equation at the interfacial 


nodes ; 




The above equation (11) is applicable only if the 
governing equation (7) holds for the layers on either side 
of the interfacial node. 


2. Finite-Difference Equation For The Ice-Abrasion 
Shield Interface 

If i is the interfacial boundary point at the ice and 
abrasion shield interface, as shown in Figure 4b, it becomes 
necessary to consider which of the governing equations 
(8a,b,c) are applicable along with equation (7) for the abra- 
sion shield and the boundary conditions (9a, b). Initially, 
the interface temperature is less than 0^-A6 and, therefore, 
equation (8a) is used. For j = 5 at node i, equation (7) 
becomes (with q''^ = 0) , 



(12a) 



For j = 6 at point i, equation (8a) reduces to 



and the boundary conditions ( 9 a,b) give: 


(12b) 




(12c) 


(12d) 


n n+1 n n+1 

The temperatures 65,1.1 and 

are fictitious and have to be eliminated between the equa- 
tions (12a,b,c,d) to give: 



( 13 a) 


n 

If the temperature at node i, 0 ^ . , is between 0 ,c"A 0 

D f r r 

and 0£+A0, equation (8c) is used along with equation (12a, c,d) 
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n n+1 

and, as before, the fictitious temperatures Sc . 0 . , 

.Of It j. 0 ^ Xt X 

n n+1 

0, . , and 0- are eliminated to give: 
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If the temperature at the interface i, 0, ., is 

b , 1 

greater than 0^+A0, equation (8b) is applied along with 
equations (12a, c,d) . Again fictitious temperatures are 
eliminated to yield the applicable equation: 



4“ 


Xu, AX5^ e?^.‘ 
(i _ C_LC^f) 
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3. Finite-Difference Equation For The Substrate- 
Ambient Boundary 

The interface under consideration is shown in Figure 
. At this grid point i equals 1. The point i equals 0 


4c 
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is a fictitious point as it falls out of the composite 
body boundary. At i equals 1, the governing equation (7) 
is applicable along with the boundary condition (5a'). 

The finite-difference representation of equation (5a') 
using equation (6a) gives, for j = i = 1, 




Similarly, the governing equation (7) reduces to 


(14a) 




e 




(14b) 


n 


Eliminating the fictitious temperatures 0^^ q, and 


n+1 

0^ Q between equations (14a,b) yields 


e 


>^+l 

iA • 


■0 


Vk 
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4. Finite-Difference Equation For The Ice-Ambient 
Boundary 

The interfacial point, as shown in Figure 4d, is the 
grid point m. The point m+1 falls outside the composite 
body boundary and is therefore a fictitious point which has 
to be eliminated. The governing equation (8a) for j = 6 
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when applied to point m is of the form: 

/ .2.^ - 2. Y>4-I 

-c.. 


= -ec.«-, 


(16a) 


As before, the boundary condition (5b'), after using 
the finite-difference representation (6a), becomes; 


^jCrOC: K.] - 


Elimination of the fictitious temperatures 6, and 

n+l 

®6 m+1 equations (16a,b) gives the finite-difference 

equation for the point m as: 






J -^5 


(16c) 


In deriving equation (16c) , it has been assumed that at 
the upper boundary, point m, ice is present and the governing 
equation (8a) is applicable. However, if instead of ice a 
"mushy" phase exists at this boundary, the governing equa- 
tion (8c) is applied at point m to yield; 


(17a) 


, *V\ A» I f“ 
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Equation (5b') after rearranging and using equation (6a) 
gives : 


A 




1 ,L 


' Vt+l v% \ 

z 


1 (17b) 


n+1 


Elimination of the fictitious temperatures 0 


n 


6 ,m+l 


and 


®6 m+1 equation (17a, b) gives the finite-difference 

equation for the point m as: 
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Finally, instead of ice or a "mushy" phase, water is 
present at the upper boundary, the governing equation (8b) 
is applied at point m to give: 




z, - A 


.»4-l 
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Applying the boundary condition (5b') in its finite- 


difference form: 
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n rv''' 

r."\ ol 

\L fi.H + 

n 

and eliminating the fictitious temperatures 6^ and 

n+1 

Sg between equations (18a, b) gives the finite-difference 
equation for the point m as: 


"!■ f^'~ -?n2,L 0 — Z. A.x^0, 

L\ ^r- ^ 

Ua J Va) 


(18c) 


5. Finite-Difference Equation For a Point Source Posi- 
tioned At An Interface Between Two Layers 
If the heater is not of finite thickness, but is 
instead a point source at node i between slabs j and j+1, 
the finite-difference equation is obtained by using the 

th 

governing equation (7) at point i for both the j and j+1^ 
slab along with the boundary conditions: 


• = 0 ■ For all n (19a) 

=- — For all n (19b) 

-r>M. 


The above equation in finite-difference form becomes: 
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Equation (7) for the j and j+1 layers, when applied at 
node i without the internal heat generation term gives: 


0 
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n n+1 n n+1 

The temperatures 0. 0. .,,, 0.,, . , and 0 . , , . , 

3,1+1' 3,1+1' 3 + 1, 1-1 3 + 1, 1-1 

are fictitious temperatures and are eliminated between the 

equations (19a,c,d,e), giving the finite-difference equation 

applicable for a point heat source as: 



( 20 ) 
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IV. METHOD OF SOLUTION AND 
PROGRAM ALGORITHM 


A. METHOD OF SOLUTION 

The numerical formulation of the de-icer pad using the 
finite difference scheme results in a set of linear equa- 
tions which can be put in tridiagonal matrix form. The 
equations are shown below and are readily solved using the 
Gaussian Elimination method suggested in Reference (31) . 


^ 1®1 ^ 1®2 


= ‘^l 

^ 2®1 ^ 2®2 

^^ 2^3 

= ^2 

^ 3®2 

^ 3^3 ^ 3 % 

= ^3 


a.e..j^ + b.e. + 0.0.^^ 

= d. 
1 


( 21 ) 


N N-1 N N N 

The complete algorithm for the solution of the tridiagonal 
system is 

= Gi - , i =■ N-l,N-2, ,1 


where BE^ = b^^ and G^^ = d^^/BEj^ 


BE^ = b^ - (a^c^_^)/BE^_^, i = 2,3, 


Gi = (d^ - a^G^_j^)/BE^, i = 2,3, 


( 22 ) 


,N 
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In the present study, non-dimensional temperatures, 0^, 
at the time step, n+1, are calculated using the above 
algorithm. The constants d^^, d 2 ,.-»»,djg are initially 
calculated using the temperatures, 6^^, at the previous time 
step, n. 

B. NUMERICAL PROGRAM ALGORITHM 

The flow diagram of the main program for the de-icer 
pad is shown in Figure 5. The computer program is listed 
in Appendix I . 
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V. NUMERICAL RESULTS 

In order to test the resulting program and algorithm, 
an unsteady state heat conduction problem involving a two 
layered gypsum board-steel composite body was solved applying 
the procedure mentioned earlier in Chapters III and IV. The 
results are shown in Figure 6 and are in excellent agreement 
with the analytical solution given in Reference (31) . The 
above problem applied convective boundary conditions at 
both boundaries; however, the computer model in this study 
has been designed to also handle constant temperature or 
mixed boundary conditions. 

Initially, the de-icer pad problem was solved assuming 
no phase change in the ice layer and uniform heat input for 
both the point heat source and the finite thickness heaters. 
For the de-icer pad configuration shown in Figure 7, the 
temperature rise at the ice-abrasion shield interface in 
order to raise the interface temperature to 0° C is compared 
with the results obtained numerically by Stallabrass (4) and 
analytically by Campbell (3). As shown in Figure 7, excel- 
lent agreement is achieved between the two numerical methods 
and the analytical solution, up to 5 seconds. Beyond 5 
seconds, the numerical solution of Stallabrass (4) gives a 
slightly optimistic temperature rise, while the present 



method shows continued excellent agreement with the analy- 
tical solution. 

To determine the effect of the various variables on the 

de-icing time, a number of cases were investigated and, 

whenever possible, compared with the results of Stallabrass 

(4) . In Figure 7, as well as in all of the following cases, 

2 

values of h^ = 1 Btu/ft hr°F and h 2 = «> were used. 

A. EFFECT OF POWER DENSITY 

The de-icing time in the present study is assumed to be 
the time interval beginning when power is applied to the 
heater and extending up until the ice-abrasion shield inter- 
face just reaches a temperature of 32 ®F. Figure 8 shows 
the results of the effect of variation in the power density 
on the de-icing time for various ambient temperatures. Good 
agreement is obtained with the results of Stallabrass (4) . 

The slight variation occurs because of the more optimistic 
nature of Stallabrass' results as indicated in Figure 7. 

Figure 8 does demonstrate what actual tests have indicated, 
namely that the total energy required in order to shed the 
ice increases with a decrease in ambient temperature and 

with a decrease in the power density. Hence, lower power 

2 

densities on the order of 15 or 20 watts/in should not be 

used for low ambient conditions because of the long de-icing 

times that are needed. Tests have indicated that a power 

2 

density on the order of 25 watts/in is the practical minimum. 
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B. EFFECT OF INSULATION THICKNESS RATIO AND INSULATION 
MATERIAL 

An insulation thickness of 0.010 inches has been 
assumed for the outer epoxy/glass insulation for the pre- 
sent study. Since it is desired that a maximum amount of 
energy released from the heater be directed toward the ice 
layer, the outer insulation should be thinner than the 
inner insulation. Figure 9 shows the effect of varying 
the insulation thickness ratio, inner insulation thickness/ 
outer insulation thickness, on the de-icing time. As the 
ratio is increased from 1 through 5, the de-icing time 
decreases appreciably. For an ambient temperature of -25 °C, 
a reduction of 20% in the de-icing time occurs as the ratio 
is increased from 2 to 5. However, further increasing of 
the insulation thickness ratio only affects the de-icing 
time for lower ambient conditions. As shown in Figure 9 
for an ambient temperature of -25 °C, the de-icing time 
does not change for an increase in thickness ratio from 5 
to 10, although it does for lower ambient temperatures. 

A decrease in the de-icing time may also be achieved 
by changing the inner insulation material. Figure 10 shows 
the effect on the de-icing time when the epoxy/glass insula- 
tion of 0.22 Btu/hr f t °F thermal conductivity is replaced 
by the same thickness of polytetraf luoroethylene (KEL-F) 
having 0.04 Btu/hr ft °F thermal conductivity. For an ambi- 
ent temperature of -25 °C, the de-icing time is reduced by 
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approximately 36%. Hence it is advantageous to use an 
insulation material of very low thermal conductivity if it 
is also a good electrical insulator. 

C. EFFECT OF SUBSTRATE MATERIAL 

Figure 11 shows the effect on the de-icing time for 
various types of substrate materials. It might be expected 
that if the aluminum alloy is replaced by stainless steel, 
which has a lower thermal conductivity and thermal diffusiv- 
ity, the de-icing time would be reduced; however, just the 
opposite occurs. This result is attributed to the higher 
thermal capacity per unit volume of the stainless steel. 
Figure 11 also shows that if an insulation layer of 0.087 
inch epoxy/glass is added to the aluminum alloy substrate, 
thereby reducing the overall thermal conductivity of the 
substrate, a decrease in the de-icing time is observed. 

D. EFFECT OF VARIABLE HEAT INPUT 

In the previous sections, the heat input was assumed 
uniform and constant. In this section, a step-wise heat 
input is applied and the temperature response at all the 
interfacial nodes is determined. Figures 12 a and b show 
the temperature response for a variable point heat source 
(3 seconds on, 1 second off) as compared to a constant 
point heat source for an ambient temperature of 5 °F (-15 ®C) . 
As expected, the temperatures at the various interfacial 
nodes drop as the heat is switched off, and begin to rise 
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again as the heat is switched on. Figures 13 a and b show 
the corresponding results for a 5 seconds on, 1 second off, 
finite thickness heater for an ambient temperature of -22 °F 
(-30 °C) . Thus, if heat is applied such that the heater 
is on for a period until the ice-abrasion shield interface 
reaches a temperature of 32 °F, and then is switched off 
until another layer of ice forms on the blade surface, the 
total energy usage for de-icing would be reduced. The 
computer program in this study has used an arbitrary periodic 
step-wise heat input, but it can be modified to apply to any 
other type of variable heat input. 

E, EFFECT OF PHASE CHANGE 

In the earlier part of this study, as in the work of 
Stallabrass (4) , the phase change in the ice layer has been, 
neglected. Hence, optimistic de-icing times were obtained 
because the latent heat needed to melt the ice has been 
ignored. To rectify this, a method used by Bonacina et al. 
(26) is applied to account for the phase change in the ice 
layer. In order to test the accuracy of this method, the 
one-dimensional water-ice solidification problem was solved. 

As mentioned in Chapter II, the selection of the phase 
change temperature interval, 2 AT,. depends upon the physical 
problem and, in combination with the number of nodes in the 
ice layer, affects the accuracy of the solution. Figure 14 
shows that a phase change temperature interval of 2 °F 



40 


provides good agreement with the results of Bonacina et al. 
(26) . Figure 14 also indicates that increasing the number 
of nodes from 81 to 126 for the temperature interval of 
2 °F does not affect the results. For the ice layer in the 
de-icer pad problem, various numbers of nodes were used to 
determine what effect, if any, they would have on the accu- 
racy of the solution. It was observed that increasing the 
number of nodes above 60 did not change the solution and, 
therefore, for each of the solutions, 60 nodes were used 
in the ice layer. 

As in the earlier sections, the de-icing time is assumed 
to be the time at which the interfacial temperature at the 
ice-abrasion shield interface reaches 32 °F. Figure 15a is 
a repetition of Figure 7, but with the phase change being 
considered, and shows that the de-icing time is increased. 
This is expected since the latent heat of the ice is now 
accounted for, whereas in the previous cases the ice had 
been treated as a single phase. Figure 15b parallels Figure 
9 and Figure 15 c parallels Figure 11 in illustrating the 
same increasing time result. Thus, the results shown in 
Figures 15 a, b and c should therefore be more realistic. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 

The one-dimensional computer simulation model developed 
in this study for the de-icer pad configuration accurately 
predicts the temperature profiles for any type of boundary 
conditions or thermal heat sources. The results agree well 
with previous numerical calculations done by Stallabrass (4) 
for cases when the phase change is not considered. The 
method of Bonacina et al. (26) to describe the phase change 
was incorporated into the model and adequately predicts the 
thermal history of the de-icer pad when the latent heat 
effect of the ice is taken into account. 

The next study should concentrate on developing a two- 
dimensional de-icer pad model in order to investigate the 
effects of blade shape, heater gap-width and heater geometry 
on the de-icing time. 
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TABLE 1 

THERMAL PROPERTIES OF SELECTED MATERIALS 


Material 

Thermal 

Spec. 

Dif fusivity 

Density 


Conductivity 

Heat 


9 

lb. 



Btu 

cal 

Btu/lb.°F 

ft^ 

cm 

g 


hr. ft.°F 

sec. cm.°C 

Cal/g’C 

hr. 

sec. 

ft^ 

cm^ 

Aluminum (soft) 

124 

0.51 

0.23 

3.20 

0.820 

169 

2.71 

Aluminum Alloy, 

66. 5 

0.275 

0.23 

1.65 

0.427 

175 

2.80 

75ST6 








Nichrome 80/20 

7.6 

0.031 

0.107 

0.138 

0.035 

515 

8.25 

Stainless Steel 

8.7 

0.036 

0.118 

0.150 

0.0385 

495 

7.93 

304 








Epoxy /Glass 

0.22 

9 X 10“‘^ 

0.23 

0.0087 

2.2 x 10“^ 

110 

1.82 

filled 32% 








Polytrif luoro- 

0.04 

1.7 X lo"'^ 

0.216 

0.0014 

0.38x10"^ 

130 

2.08 

chloroethylene 

(KEL-F) 








Water (0°C) 

0.320 

1.32xl0~^ 

0.997 

0.0051 

0.00132 

62.4 

1.0 

Ice (pure) 0°C 

1.293 

5. 35xl0“^ 

0.5057 

0.0445 

0.0115 

57.2 

0.9168 

-10°C 

1.356 

5.61xl0"^ 

0.5038 

0.0469 

0.0121 

57.3 

0.9182 

-20°C 

1.416 

5.86xlo"^ 

0.5020 

0.0492 

0.0127 

57.4 

0.9196 


Latent Heat of Fusion of Ice = 143.4 Btu/lb. = 79.75 Cal/g. 
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1 METAL SUBSTRATE 2 INNER INSULATION 
3 HEATER ELEMENT 4 OUTER INSULATION 
5 ABRASION SHIELD 6 ICE 


Figure 1. One Dimensional De-icing Pad Model 
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Figure 2. Section Of Wing Fitted With Pneumatic Boot 
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Figure 3. One Dimensional Finite-Difference De-icer 


Pad Model 
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Figure 4. Finite Difference At Selected Grid Points 
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Figure 5. Flow Diagram Of Main Program 
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Figure 6- Comparsion Of Finite-Difference And Analytical Results ror Composite Slab 
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Figure 7. Comparison Of Finite-Difference And Analytical Results For De-icer Pad 
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Figure 8. Effect Of Power Density On De-icer Performance 
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Figure 13b. Interfacial Temperature Response To Constant 
And Variable Finite Thickness Heat Source 
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APPENDIX I 


(XIMPLETE PKX3^ LISTING 

NUMERICAL EMULATION OF ONE -DIMENSIONAL HEAT TRANSFER 
IN COMPOSITE BODIES WITH PHASE CHANGE. 

THIS PROGRAM CAN CALCULATE THE TEMPERA fURE PROFILE IN 
A COMPOSITE SLAB WHICH HAS CGNUEC T lUE f CONS (ANT TEMP- 
ERATURE OR MIXED BOUNDARY CONDITIONS. 

THE PROGRAM CAN ALSO BE USED FOR COMPOSITE BODY PROBLEMS 
WITH CONSTANT OR VARIABLE HEAT SOURCES. 


A . B » C y D 

AK 

ALP 

AKL 

ALPL 

CA 

CPSyCPL 

DTAUI 

DTAUM 

DTAUF 

DX 

EL 

HEAD 

HI 

H2 

HLAM 

IBCl 

IBCl 

IBC2 

IBC2 

ICOUNT 

I FREQ 

IG 


= TRIDIAGONAL MATRIX CONSTANTS. 

= THERMAL CONDUCTIVITIES OF LAYERS. 

= THERMAL DIFFUSIVI I lES OF LAYERS. 

= THERMAL CONDUCTIVITY OF WATER. 

= THERMAL DIFFUSIVITY OF WA I ER . 

= CONSTANT. 

= SPECIFIC HEAT OF ICE AND WATER PER UNIT VOLUME. 
= INITIAL TIME STEP 
= INTERMEDIATE TIME STEP 
= FINAL TIME STEP 
= SPACING BETWEEN NODES. 

LENGTH OF EACH LAYER. 

HEADINGS. 

= HEAT TRANSFER COEFF. AT LOWER BOUNDARY. 

= HEAT TRANSFER COEFF. AT UPPER BOUNDARY. 

= LATENT HEAT OF ICE. 

1» IMPLIES TEMPERATURE IS CONSTANT AT X=0. 

= 2» IMPLIES CONVECTIVE HEAT TRANSFER AT X-0. 

= ly IMPLIES TEMPERATURE IS CONSTANT AT 
= 2f IMPLIES CONVECTIVE HEAT TRANSFER AT X=l. 

= COUNTER ON TIME SIEP. 

= NUMBER OF TIME STEPS BETWEEN SUCCESSIVE 

PRINTING OF THE TEMPERATURE PROFILE. 

= ly IMPLIES PHASE CHANGE IN ICE LAYER IS NOT 

CONSIDERED. 
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IG 

= 2f 

IMPLIES 

IH 

= If 

IMPLIES 

IH 

= 2 r 

IMPLIES 

IH 

= 3f 

IMPLIES 

IHQ 

= If 

IMPLIES 

IHQ 

= 2f 

IMPLIES 


PHASE CHANGE IN ICE LAYER IS 

NO HEAT SOURCE. 

POINT HEAT SOURCL.. 

HEAT GENERATION WITHIN SLAB. 
CONSTANT HEAT SOURCE. 
VARIABLE HEAT SOURCE, 


CONSIDERED 


IJ 

INTIME 

IMTIME 

M 


SLAB WITHIN WHICH HEAT GENERATION OCCURS. 
NUMBER OF TIME STEPS FOR WHICH INITIAL TIME 

STEP IS USED. 

NUMBER OF TIME STEPS FOR WHICH INTERMEDIATE 

TIME STEP IS USED. 

NUMBER OF NODES IN SLAB. 


MM = INTERFACE NODE NUMBERS. 

N = NUMBER OF LAYERS IN SLAB. 


N1 

N2 

NO 

NOl 

N02 

NODE 

Q 

cr2 

OHEAI 

QV 

T 

TDIFF 

IE 

TGI 

TG2 

TIN 

TLEN 

TMAX 

rOFF 


= LOWER SLAB NUMBER FOR POINT MEAT SOURCE. 

= UPPER SLAB NUMBER FOR POINT i lEA X SOURCE, 

= NUMBER OF NODES IN EACH LAYER. 

LOWER NODE NUMBER TOR FINITE THICKNESS HEATER. 

== UPPER NODE NUMBER FOR FINITE THICKNESS HEATER. 

= NODE AT WHICH POINT HEAT SOURCE IS APPLIED. 

POINT HEAT SOURCE WATTS/ INH<IN, 

= VOLUMETRIC HEAT SOURCE WA TTS/IN)iaN*IN . 

= FUNCTION PROGRAM FOR VARIABLE HEAT INPUT. 

STEP INPUT FOR VARIABLE HEAT SOURCE. 

= NON-DIMENSIONAL TEMPERATURE. 

= HALF PHASE CHANGE lEMI'ERATURE INTERVAL. 

= TEMPERATURE. 

= AMBIENT TEMPERATURE AT LOWER BOUNDARY OF SLAB. 
= AMBIENT TEMPERATURE AT UPPER BOUNDARY OF SLAB. 
= INITIAL TEMPERATURE IN SLAB. 

= TOTAL LENGTH OF SLAB. 

= ICE-ABRASION SHIELD INTERFACE TEMPERATURE. 

= OFF TIME OF STEP HEAT INPUT. 


ION = ON TIME OK STEP HEAT INPUT. 

rPHAS ^ TEMPERATURE AT WHICH PHASE CHANGE OCCURS. 



67 


TR 

TREF 

TRIDAG 

TXO 

TXl 


DATA IN/5/»I0/6/ 

DIMENSION HEAD(40r80) 

DIMENSION A( 100) »B( 100) rC< 100) r D ( 100 ) » T ( 200 r 100 ) yTRdOO) i-TE( 100) 
DIMENSION ALP(6) » CA ( 6 ) » NO ( 6 ) r MM ( 6 ) » DX ( 6 ) » AK ( 6 ) y EL ( 6 ) 

INPUT DATA 

DO 300 1=1 r 18 

READ ( IN y 100) (HEAD( I y J) y J=1 ySO) 

300 CONTINUE 

READdNylODNyMy FLEN 

DO 10 I=lylOO 
A(I)=0, 

B ( I ) =0 . 

C(I)=0. 

D(I)=Oy 
TEC I )=0. 

TRC I )=0. 

DO 11 J=ly200 
T( Jyl)=0. 

11 CONTINUE 
10 CONTINUE 
DO 12 K=lyN 
ALP(K)=0. 

AKCK)=0. 

£LCK)=0. 

DX(K)=0. 

CACK)=Oy 


= NON-DIMENSIONAL TEMPERATURE AT PREVIOUS 
= REFERENCE TEMPERATURE y 

= SUBROUTINE TO SOLVE TRIDIAGi’UL MATRIX* 

= CONSTANT TEMPERATURE AT LOWER SLA*B BOUNDARY* 
= CONSTANT FEMPERATURE AT UPPER SLAB BOUNDARY* 
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N0(K)=0 

MM(K)=0 

12 CONTINUE 

INPUT DATA 
DO 301 I=1V,22 

301 READ ( IN f 100) <HEAD( X , J) » J~1 ,80) 

DO 13 K=1,N 

READ ( IN, 102) NO(K) ,EL(K) , AK ( K ) , ALP ( K ) 

13 CONTINUE 

DO 302 1=23, 25 

302 REAtK IN, 100) (HEAIK I , J) , J=1 ,00) 

READ; IN, 103) III, NODE, N1 ,N2 
read; IN, 10^1 ) I J,N01 ,N02, IHO 
read; IN, 105)0, TON, TOFF, 00 

DO 303 1=26,20 

303 read; IN, 100) ;head; 1 , J) , j=i ,eo) 
read; IN, 106) IBCl , [BC2 

read; in, 116 ) rxo, rxi 
read; IN, 107 ) TGI , HI . IG2,H2 
read; IN, 108) riN, irlf 
do 304 1=29.31 

304 READ ; IN , 100 ) ( HEAD ; I , J ) , J=1 , 80 ) 

RE AD ; I N , 1 09 ) I G , I II. AN , ALI ‘L , AKL , CPL , CPS 
read; IN, 110)TPHAS, TDIFF 
DO 305 1=32,34 

305 read; IN, 100) ;head;i, J) , j=i,bo) 

read; in, 1 1 1 ) DTAUI , INSTEP , DTAUN, INSTEP 
READ ; I N . 1 1 2 ) DTAUF , IFREQ , IMAX 
DO 315 1=35,40 

315 read; IN, 100) ;head;i,j) ,j=l,oo) 

PRINT THE DATA 


DO 306 1=1,18 



306 WRITE( I0f900) (HEAIK I» J) . J=1 »80) 

URITE( I0»200)N»f1» FLEN 

DO 307 1=19 » 22 

307 WRITE( IQ»900) (HEALK I » J) , J=1 f30) 

DO 14 K=1»N 

URITE(I0»102)N0<K) » EL ( K ) » AK ( K ) »ALP(K.) 
14 CONTINUE 

DO 300 1=23 » 25 

308 UR1TE( 10 » 900; <HEAD( I » J) » J=1 ,80) 
IFdH.EQ.l) GO TO 60 
IF<IH.£Q.2) 00 TO 61 

WRITE( 10,201 ) IJrNUl ,N02 
GO TO 62 

60 WRITE ( 10,202) 

GO TO 62 

61 WRITE( I0,203)N0DE,N1 ,N2 

62 [ F (I HQ . EQ . 1 ) WR I TE < 1 0 , 204 ) Q 

IF (IHQ.EQ. 2) WRIT E( 10,224) TON, TOFF, QV 
DO 309 1=26,28 

309 WRITE ( 10,900) (HEAD(I,J) , J=l,80) 

IF< IDCl ,EQ.2) GO lO 63 

WRITE ( 10,205) rxO 
GO TO 64 

63 WRITE( I0,206)TGl ,H1 

64 IF( IBC2.£a.2) GO TO 65 
WRIT£(I0,207)TX1 

GO TO 66 

65 WR I TE< 10,208) TG2,H2 

66 WR I T E < 1 0 , 209 ) I I N , T RFJ- 

WR I TE( 10,900) CHEAD(26, J) , J=1 ,80) 

IF( IG.EQ. 1 )G0 TO 67 
WRITE( 10,210) 

WRITE ( 10,900) (HEAD<29, J) , J=l,80) 
WRITE( 10,900) (HEAD(28, J) , J=l,80) 
WRITE( 10,900) (HEAD(30,J) ,J=1,80) 
WRITE( 10,900) (HEAD(31, J) , J=l,80) 


70 


URITE( IQ»211)HLAM» AKL» ALPL»CPL»CPS 
URITE( I0»212)TPHAS,TDIFF 
GO ro 68 

67 WRITE ( 10,213) 

68 WRITE (IQ, 900) (H£AD(26, J) , J=l,80) 

DO 310 1=32,34 

310 WRITE( 10,900) (HEAD(I,J) ,J=1,80) 

WRITE( ia,218)DTALII , INSTEP , DTAUM , IfiSTEP 
WRITE( I0,214)DTAUF, IFREa,TMAX 
WR I TE( 10,900) (HEAD(1,J) ,J=1,80) 

WRIT£( 10,217) 

WR I TE ( 1 0 , 900 ) ( HE AD < 1 6 , J ) , J= 1 , BO ) 

URITE< 10,900) ( l-IEAD( 1 ,J) ,J=1 ,80) 

INITIAL CONDI r IONS 

FD IST=TIN/TREF 
G1TIME=TG2/TR£F 
G0TIM£= I Gl/ IREF 

xitime=txi/tri:f 

XOTIME=TXO/ [RLE 
rLEN=TLEN/12. 

DO 30 J=1,N 
30 EL( J)=EL( J)/12. 

T0N=TGN/3600. 

TOFF=rOFF/3600. 

DTA0I=DTA0I/3600 . 

DTA0h=DTA0M/3600 . 

DTA0F=DrA0F/3600. 

INTERFACE NODE NO.iDERING 

flM< 1 )=N0( 1 ) 

DO L6 L=2,N 
MM<L..)=Mh(L-l)iNO(L)-l 
16 CONTINUE 


■?! 
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PRINT THE INITIAL lEMPERATURE PROFILE 

ICOUNT=l 
no 19 J=l.M 
FLOAT J=J 

19 T< ICGUNTr J)^=FDIGT 
TAU=0. 

URITEC tOf 21S) TAIJ 
no 31 I=1»M 

31 TE ( 1 ) = r < ICOUNT * I ) tTREF 
WRITE ( 10^216) < IE( I ) r I---^! 

88 LCOUNT=ICOUNT 
ICOUNT^LCOUN r f 1 

IF ( ICOUNT. LE. INSTEP ) nr AlJ-LiTAU I 
I r S T E P = I N S r E F-' f I M S r Fi. P 

IF ( ( ICOUN 1 . LE . I TS T El ' , ANH , ( ICOUNT . GT . INSTEP ) ) U I AU^^DTAUM 

IF ( ICOUNT . GT , I TSTEP ) nTAU==OTAUF 

MCL=200 

IF ( ICOUNT. GT ,hCL) GO lU 999 
TAU=TAUfn lAU 

CALCULATION OF CONSIANTS 
no 15 1 = 1, N 

nx ( I ) =£L ( [ ) / ( TLEN* ( NO ( I ) -1 ) ) 

15 CA( I ) = ( ALP( I ) ^KThU) / <nx( I )>KnX( l)*rLEN* TLEN) 

CALCULATION OF THE IRiniAGNOL CUNSTANIS 

IFtN.NE.1) go to 32 
no 29 J=1,M 
A< J)=l . 

B( J)=-2.*( 1 . F-1 ./CA( 1) ) 

29 C(J)=1. 
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GO TO 34 
32 LIST=1 
NEW=N-1 
DO 17 K=lfNEW 
NN=MM(K) 

B(NN)=-1.)K< <1. fl./CA(K) ) + ( (AK(Kfl)*DX(K) )/ 
1 <AK<K)*DX<KH ) ) )*( 1 ,M ./CA(Kfl) > ) 

C(NN) = (AK<Kfl )*DX(K) )/(AK(K)*DX<KH ) ) 
A(NN)=1 . 

NNEW=NN-1 

DO 18 J=LIST»NNEW 

A( J)=l. 

B( J)=-2.#< 1 ♦ H ./CA(K) ) 

C< J) = l . 

18 CONTINUE 
L1ST=NN+1 
17 CONTINUE 


ACCOUNTING FOR PHASE CHANGE IN ICE LAYER 


34 IF(IG.EQ.l) GO TO -10 
THETA=TDIFF/TREF 

CSTAR= < ( CPSfCPL ) /2 . ) THLAM/ ( 2 . *TRCF*THETA ) 

CONST 1 = ( TLEN*TLEN*DX ( N ) )KDX ( N ) ) / ( hI...P ( N ) KDTAU ) 

C0NST2= ( TLEN* TLEN*DX ( N ) icDX ( N ) ) / ( ALPL^DTAU ) 

C0NST3=2 ♦ *TLEN*TLEN)tcDX ( N ) »DX ( N ) *C3TAR/DTAU 

THETAl=<TPHAS/TR£F)+rHETA 

THETA2= ( TPHAS/TREF ) -THETA 

TKK 1 = < AKL-AK < N ) ) / ( 2 . JtcTHETA ) 


ACCOUNTING FOR PHASE CHANGE IN THE ICE LAYER 


40 IF(N.EQ.i) GO TO 35 
IF( IG.EQ.l )G0 TO 42 
LUST=LIST-1 



CALCULATION OF THERMAL CONDUCTIVITIES AND TEMPERATURES 
AT INTERMEDIATE TIME STEPS FOR •MUSHY* REGION. 

DO 20 J=LUSTrM 
IF( J.EQ.LUST) GO TO 401 
IF(T(LCOUNr, J) .GT. THETAl) GO TO 43 
IF(T<LCOUNT» J) .LT. rHETA2) GO TO 44 

TEMPERATURE AND THERMAL CONDUCTIVITY AT NODE I TOR AN 
INTERMEDIATE TIME STEP. 

401 TKMIN1=AK(N)+TKK1*< < ( T < LCOUN T . J ) TT ( LCOUNT . J- 1 ) )/2. ) 
l-( (TPHAS/TREF) -THETA) ) 

IFCTKMINl .LT.AKL) TKMIN1=AKL 
IF(TKMIN1 .GT.AK(N) ) TKM INI =AK ( N ) 

T(LCOUNT»M + l)=iGl TIME 
T(LC0UNT»MT2)=G1TIME 

PLUSKl=AK<N)+TKKl:»t( ( < T ( LCOUN T r J ) +T ( LCOUNT » J+ 1 ^ )/2. ) 

1- ( (TPHAS/TREF) -THETA) ) 

IF < PLUSKl . LT . AKL ) PLU3K 1 :=AKL 
IF(PLUSK1 .GT.AK(N) ) PLUSK 1-AK ( N ) 

TEMP1 = T ( LCOUNT » J ) f ( DTAU/ ( 2 . *CSTAR* TLEN* TLEN*DX ( N ) :):DX ( N ) ) ) 
!:♦(( < PLUSKl* ( T (LCOUNT r J Tl ) -T(LCOUNTrJ) ) - TKMINI * ( T ( LCOUN C y J ) ) 

2- T(LC0UNTt J-1 ) ) ) 

TEMPERATURE AND THERMAL CONDUCTIVITY AT NODE IM FOR 
INTERMEDIATE TIME STEP. 

TKMIN2=AK(N) fTKKl*( ( ( T ( LCOUNT y J ) IT ( LCOUN T r J+ 1 ) )/2. ) 
l-( (TPHAS/TREF) -THETA) ) 

IF ( TKMIN2.lt .AKL) TKMIN2=AKL 
IF(TKMIN2.GT.AK(N) ) TKMIN2=AK ( N ) 

PLUSK2=AK(N)+TKK1*( ( ( T ( LCOUNT y JT2 ) ET ( LCOUNT y JIl ) )/2. ) 

1 - ( ( TPHAS/TREF ) -THETA ) ) 

IF ( PLUSK2 . LT . AKL ) PLUSK=AKL 

IF ( PLUSK2 . GT . AK ( N ) ) PLUSK2=AK ( N ) 
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TEMP2=T < LCOUNT » J+ 1 ) + < DTAU/ ( 2. *CSTAF<)KTLEN)tc rLEN*DX ( N ) >|cDX < N ) ) ) 

1*( <PLUSK2*(T(LC0UNT» Jf2)-T(LC0UNT» Jf 1) )~TKMIN2»(T<LC0UNT» JH ) ) 
2-T(LC0UNT» J) ) ) 

TEMPERATURE AND THERMAL CONDUCTIVITY AT NODE I -1 FOR 
INTERMEDIATE TIME STEP, 

rKMIN3=AK<N) fTKKl*( ( (T(LCUUNr» J-2)+T<LC0UNTr J--1 ) )/2, ) 
l-( (TPHAS/TR£F)-THETA) ) 

IF(TKMIN3.Lr,AKL)TKMIN3=AKL 
IF ( TKM IN3 . GT . AK ( N ) ) TKM I N3=AK ( N ) 

PLUSK3=AK(N)+TKK1*( ( ( f < LCOUNT » J ) fl ( LCOUNT , J-1 ) )/2, ) 

1- ( (TPHAS/TREF) -THETA) ) 

IF < PLUSK3 . LT , AKL ) PLUSK3=AKL 
IF<PLUSK3.GT.AK(N) ) PLUSK3=AK ( N ) 

TEMP3=T < LCOUN F , J-1 ) f ( DTAU/ ( 2 , )KCGTAR*TLEN)»tTLEN1!DX ( N ) *DX ( N ) ) ) 

!»( (PLUSK3)»c< KLCOUNF , J)-F(LCOUNT, J- 1 ; )- IKM1N3* ( I ( LCOUNT , J-1 ) ) 

2- T ( LCOUNT tJ--2) ) ) 

CALCULATION OF THERMAL CONDUCFIVITY FOR “MUSHY" REGION, 

TKMIN = AK(N) + FKK1;K( ( < lEMPl f l EMP3)/2, ) 
l-( (TPHAS/TREF) -THETA) ) 

IF(TKMIN,LT.AKL)TKMIN=AKL 
IF ( TKMIN , GT . AK ( N ) ) TKMIN=AK ( N ) 

PLUSK=AK(N)+TKKl;((( ( ( FEMP1HEMP2 ) /2 . ) 
i-( (TPHAS/TREF) -THETA) ) 

IE ( PLUSK , L T . AKL ) PLUSK^=AKL 
I F ( PLUSK . GT . AK ( N ) ) PLUSK-AK ( N ) 


THERMAL CONDUCTIVITY AT THE ICE-AMBIENT BOUNDARY 
IF(J.NE.M)GO TO 403 

PKB=AK(N) ITKK1)*:(TEMP1 -( ( TPHAS/TREF ) -THE TA ) ) 

IF(PKB,LT,AKL)PIB=AKL 

IF ( PKB . GT . AK ( N ) ) PKB=AK ( N ) 





IF( J.EQ.fDGO TO 45 


THEFiflAL CONDUCriVITY OF -MUSHY" REGION AT .1 Ci;;. -ABRASION 

SHIELD INTERrAOE.* 

403 IF( J.NE.LUST)GO TO 402 

PKINC=AK < N ) f TKKDlc ( TEMPI - ( ( TPHAS/TREF ) -THETA ) ) 

IF ( PK INC . LT . AKL ) PK I NC=AKL 
IF(PKINC.GT.AK(N) )PKINOAK(N) 

TKINC2=TKMIN 
PKINC2=PLUSK 
GO TO 20 

APPLYING HEAT LOUATION FOR * MUSHY* REGION, 

402 A(J)=TKMIN 

B ( J ) =-l . * ( TKMIN+PLU3K+C0NST3 ) 

C(J)=PLUSK 

D < J ) =-T < LCOUNT » J- 1 ) *TKM IN-PLUSKiXT ( LCOUNT > J M. ) f T ( LCOUNT i- J ) 
l)|c( TKMINfPLUSK- C0NST3) 

GO TO 20 

APPLYING HEAT EQUATION FOR LIQUID REGION. 

43 IF(J.EQ.M)GO TO 45 
A(J) = 1. 

B( J)=-2.;«c( 1 , fC0NST2) 

C(J)=1. 

D ( J ) =-T ( LCOUNT » J f 1 ) -T ( LCOUNT » J- 1 ) +T ( LCOUNT » J ) 

1*2.*(1 .-1 .*C0NST2) 

GO TO 20 

APPLYING HEAT EQUATION FOR SOLID REGION. 

44 IF(J.EQ.M)GO TO 45 
A( J)=l. 
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B( J)=-2.*( 1 , + l ./CA<N) ) 

C( J) = l . 

D ( J ) =-T ( LCOUNT » J f 1) - T ( LCOUNi » J-1. ) +T ( LCOUNT r J ' 
- l./CA(N) ) 

20 CONTINUE 
GO TO 45 

CALCULATION OF CONSTANT FOR SINGLE SLAB. 


42 A(M)=i. 

C(M)^=1 . 

MEU=M-1 

no 33 J=LIST»MEU 


A ( J ) =•• t . 

C < 3 T = I . 

B( J) = -2.»( 1 . H ./CA(N) ) 

D ( J ) =-T ( LCOUN T r J M ) "T ( LCOUN I , J-1. ) +T ( LCOUNT » J . 
1*2. *( 1 .-1 ./CA(N) ) 

33 CONTINUE 
45 LA=2 
A < M ) 1 . 

C'. 


CALCULATIONS OF THE INTERFACE CONSTANTS 


NEW^-N-1 

DO 21 K=lrNEU 

NN=NM(K) 

n(NN)--T (LCOUNT »NNT 1 ) * (H AK ( K + 1 ) *nx ( K ) ) / ( AE ( K ) *DX ( K f 1 > ) ) 
i f < ( 1 . -1 ./CA(1\ ))♦•(< AE(K+1 >*nX(K^ >/( AK(K^*LiX (K !■ !, ; ) ) 

2*(1 . -1 ./CA(K f-1 ) ) ) * T (LCOUNT»NN) "KLCOUNTf NN-1 . 

NNEU=NN-1 
DO 22 J=LArNNEU 

D(J)=-T( LCOUNT rJ + 1 ) -KLCOUNTr J-1 ) IT ( LCOUN T . J ) 
l*(2.-2./CA(K ) ) 

22 CONTINUE 



56 D(M)=-T(LCOUNrrM-l ) - (2.*rLEN;t:H2*DX(N)/AK(N) ))KGl'nME 
1 H ( 1 . - 1 . /CA ( N ) ) f ( H2»TLEN*DX ( N ) /AK ( N ) ) ) ♦ F ( LCOLJNT r M ) 

B ( M ) =~ ( 1 . f 1 . /CA ( N ) ) - < H2* TLE;N*DX ( N ) ) /AK ( N ) 

A(M)=1, 

GO TO 49 

48 

D(MIN)=D(rtIN)-XlTIME»C<M [N) 

T( ,[eOUNT»M)-X.l 1 ln£ 

CALCULATION AT FHE LCE-ABRAS I ON SHIELD INTERFACE 

49 INC=MM(N-1) 

[F(IG.EQ.l) GO ro 50 

IF(T(LCOUNTf INC) .LF . IHLTA2)G0 FG 50 

IFCKLCOUNTr INC) .OF. FHETAI ) GO TO 51 

DUO T = DX ( N-1 ) *PK INC/ ( 2 . *TK INC2:tcDX ( N ) ^AK ( N-1 ) ) 

A( INC) = 1 . 

BaNC)=-( 1 . I I ,/CA(N-l ) ) -(DUCT*( FKINC2 + PKINC2 + C0NST3) ) 

C (I N C ) = D U C T * ( F K I N C 2 f P K I N C 2 ) 

D(1NC )=-DlJCF*(TKINC2 H-'K 1NC2 ) F ( LCOUNT r INCH )-T (LCOUNT y INC-1 ) 
1 \ ( ( 1 . -1 ./CA(N-1 ) ) I ( FKINC2TPK1NC2 -CCNST3).KDUCT);|CF(LCGUNTyINC) 
GO TO 50 
51 A(INC)=1. 

B< INC )=-l .;♦:(( 1 . l-l . /CA( N-1 ) ) H ( AKL^DX UJ -1 ) ; / ( AK UN-1 ) ;KBXFN) ) ) 
!*(!. + ( ALP ( N ) / (C A ( N ) ALPL ) ) ) ) 

C( [NC)^(AKL*DX(N 1 ) )/<AK(N-l );>!DX(N) ) 

B ( I NC ) =- ( AKL*DX ( N- 1 ) / ( AK ( N- 1 ) >|:FiX ( N ) ) ) % J ( LCOUNT i. INC+ 1 ) 

1 - F ( LCOUNT y INC- 1 ) I T ( LCOUNT y INC ) * ( ( 1 . - 1 . /CA ( N - 1 ) ) T 

2 (1 . - ( ALP ( N ) / ( ALPLiKCA < N ) ) ) ) * ( AKL-KDX ( N-1 ) / ( AK ( N-1 ) >KDX ( N ) ) ) ) 

HEAF FERN INCLUSION 

50 1F( IHO.EQ. 1 )Q1=FF 

IF ( II IN . EG . 2 ) 01 =^NHEAT ( TAU y TON y FUFF y 00 ) 

IF(IH.EO.l) GO FO 52 
IF(tH.E0.2) GO TO 53 


% 
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LA=NN+1 
21 CONTINUE 

ACCOUNTING FOK' BOUNDARY CONDITIONS 

35 IF(N.NE.i) GO TO 36 
MEW==M--1 

DO 37 I=2»hEU 

37 n ( I) = - r ( LCOUN r » I n ) -T ( lcount » i - 1 ) + r ( lcount » i ) 

l*(2.-2./CA( 1 ) ) 

36 IF- ( IBCl .EO. 1 ) GO TO 46 

D( 1 )=-T(LC0UN1 »2)-(2.*H1;KTLENI<DX( 1)/AK< 1 ; ):KG0TINE 
1-H ( 1 . -1 ,/CA( 1 ) )KH1* rLEN:4< (:iX( 1 ) /AK( I ) ) ) ♦ I ( LCOUNT i- 1 ) 

B( !)=-( 1 . H ./CA( 1 ) )-(FU*TLEN*DX(1. ) /AK( 1 ) ) 

ca)=i. 

GO TO 47 

46 D(2)=D(2)-X0TIME 
T( ICOUNT* 1 )=XOTIME 

47 IF( IBC2.E0. 1 ) GO TO 48 
IF(IG,EO,UGO rO 56 

ACCOUNTING TUIM I 'I I AGE CHANGE AT UF'PER BOUNDARY 

IF(T(LCOUNITN) .LI . niETA2)G0 TO 56 
IF(T(LCOUNITM) .Gf. THETADGO TO 57 
A(M)=TKMINfPLUSK 

B ( h ) =- ( ( PLUGK-4- TKM I N-TC0NST3 ) + (. 2 , *H2*TLEN*PLUSK*DX ( N ) /PKB ) ) 

D ( M ) =- ( TKM INTPLUSK ) *T < LCOUNT y M-- 1 ) - ( 4 . *H2*TLEN*PLUSK)XDX ( N ) * 
IGITIME ) /PKB FT ( LCOUNT y M ) * ( PLLiSKFTKM I N- CONS I 3F 
2(2. *H2*TLEN»PLUSK*DX ( N ) /PKB ) ) 

GO rO 49 

57 D ( fi ) =- r ( LCOUNT r r1- 1 ) - ( 2 . ♦ ILEN>KH2*DX ( N ) /AKL ) *G1 T INE 

1 M ( 1 .-(TLEN*TLEN»DX(N)*DX(N)/(ALPL*DTAU) ) )F(H2*TLEN*DX<N)/AKL) ) 
2»r (LCOUNT y/1) 

B ( M ) ==- 1 . - ( TLEN* TLEN^DX ( N ) ^DX ( N ) / ( ALPL*DTAU ))-(. H25f<TLEN>KDX ( N ) /AKL ) 
GO TO 49 



FINITE THICKNESS HEATER 


02=3.4121*144, ♦01/ ( TLEN*i:iX( IJ) ) 

H£AT=( (DX( IJ)*rLENJ**2)*02/(TREF*AK( IJ; ) 

FACT1=AK< IJ)*DX< IJ-1)/<AK< IJ-1 )*DXaj) ; 

FACT2=AK< IJ+1 )*DX( lJ)/< AK( IJ)*DX< IJH. ; ) 

IF< (N02-N01) .GT. DGO TO 54 
A(N01)=1. 

C(N01 )=FACT1 
A(N02) = 1 . 

C(N02)=FACT2 

B(N01 )=-l .*( ( 1 ,M ,/CA< [ J-1 ) ) fFACTl*< 1 . fl ,/CA( IJ) ) ) 

Ei(N02)=-l .*( ( 1 .H ./CA< IJ) ) fFACT2*(1. . H ./CA( IJfl ) ) ) 

54 EKNOl )=-T(LC0UNT»N01 -1) -FACTl*TU..C0UNT»N01fl ) 

l + ( ( 1 . -1 ./CA( I J-1 ) ) + (l . -1 ,/CA< IJ) )*FACTl)*r(LCn!JNTi.N01 )-FACTl*HEAT 
D(N02)=-T(LC0UNT rNQ2-l) -FACT2*T ( LCOUNT f N02 + 1 ; -i ,EAT 
IK ( 1 . -1 ,/CA( IJ) )K 1 .~1 ,/CA( I Jil ) )*FACT2)*T(LCUUNT»N02) 

IF( (N02-N01) .LE.DGO CO 52 

N0U=N01H 

NN0U=N02-1 

DO 55 IL = NOU.NN()U 

55 D( IL)=D(IL)-2.*HEAr 
GO TO 52 

POINT HEAT SOURCE 

53 D(NODE)=D(NODE) -( 2 . *TLEN*DX ( N1 ) *01*3 . 4121*144 . ) / ( TREF*AK ( N1 ) ) 

SOLOING THE TRIDIAGNOL MATRIX 

52 IF( ( IBCl .£0.2) . AND. ( IBC2.EQ.2) )G0 TO 94 
IF( (IBCl.EQ.l) .AND. (IBC2.EQ.1) )G0 TO 95 
IF( ( IBCl .E0.2) .AND. ( IBC2.EQ. 1 ) )G0 TO 96 
CALL TR I DAG ( 2 r M , A i- B » C Ki » TR ) 

TR( 1 )=XOTIME 
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GO TO 97 

95 MIN=M-1 

CALL TR I DAG ( 2 i- MIN » A . B » C » D » TR ) 

TR(1)=X0TIME 
TR<M)=X1TIME 
GO TO 97 

94 CALL TRIDAGC 1 »M»A,B»C»D»TR) 

GO TO 97 

96 M1N=M-1 

CALL TRIDAG( 1 *MIN.A»B»C»n,TR) 

TR<M)=X1TIME 

97 TMAX1=TMAX/TREF 
DO 23 J=1.M 

T< ICOUNT» J)-=rR( J) 

TE( J)=TR( J)*TREF 

23 CONTINUE 

IF( (ICOUNT/IFREQ)*IFR£O.NE.ICOUNT)GO TO 38 

M0UER2=INC 

TITAU=TAU*3600. 

WRITE<IO»215)TITAU 
WRITE<I0»216) <TE(I) »I = 1.»M) 

URITE(IOfVOO) (HEAIK 1 yj) »J=1 /BO) 

URITE( 10^900) (HEAl.i(35» J) » J=lyOO) 

UR I TE ( 1 0 » 900 ) ( HEAD (1 > J ) » J= ;i , 00 ) 

JI = 1 

NEU=N-1 

DO 24 L=1»NEU 

:i;d=35+l 

INT==MM(L) 

WRITE( lOfVOO) (HEADdD* JD) , JD=1 rOO) 
URITE<I0»216M lEd) rI = JlfINT) 

JI=INr 

24 CONTINUE 

IF(T< IC0UNTfM00ER2) .LE.TMAXl ) GO TO 08 
100 FORMAT(OOAl) 

900 FORMAT (' dOOAl) 


ro r-j rj rJ ro 
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1F13.6» 'DEG.F V/ ' THE REFERENCE TEMPERATURE ' r 20X y ' TREE =='yF13*6y 
2'DEG*F' ) 

10 FORMAT <// ' THE PHASE CHANGE IN THE ICE LAYER .1.3 CONSIDERED ' ) 

11 FORMAT (/ ' LATENT HEAT OF ICE 'y21Xy' HLAM = ' ■ i 1 3 . ' B * T . U . /LB ' / 

1' THERMAL -CONDUCTIOITY OF WATER 1 2X AKL = ' » F 1 3 . 6 f ' B . T » U ♦ /HR ♦ 
2FT.DEG.FV ' THERMAL DIFFUSIOITY OF WATER '»12X*'ALPL :^'yF13,6y 
3' FT. FT/HR,'/ ' SPECIFIC HEAT * DENSITY OF WATER '»OXi-'CPL ~'y 
4Fl3.6f 'B.T.U./CU.FT.DEG.r- V ' SPECIFIC HEAT * DENSITY OF ICE 'y 
SrlOXy 'CPS »FI3.6y ' D . T . U . /CU . FT . DEG . F ' ) 

12 FORMAT(/ ' PHASE CHANGE TEMPERATURE ' r 24X »' TPHAS ^ » F 13 . 6 y ' DEG . F ' / 

1' HALF PHASE CHANGE TEMPERATURE INTEROAL TDIFF^:: ' y FI 3 . 6 y 

2'DEG,F') 

13 FORMAT (// ' THE PHASE CHANGE IN THE ICE LAYER IS NOT CONSIDERED') 
18 FORMATC/ ' INITIAL TIME STEP ' y 26X y ' DTAU I = ' . F 1 3 . 6 y ' SECS ' /28X y 

I'FOR TIME STEPS INSTEP ^='yI3/ ' INTERMEDIATE TIME STEP'ySlXy 
2'DTAUM yF13.6y 'SECS'/20Xy 'FOR TTME STEPS IMSTEP ='yI3) 

214 FORMAT (/ ' FINAL TIME STEP ' y 2BX y ' DTAUF ' y Fi;:.: , 6 y ' SECS ' /IX y 
1' FREQUENCY OF TIME STEP/PRINT OF OUTPU T ' y 8X y ' IFREQ ='yI3/ 

21Xy' MAX. TEMP. AT ICE-ABRASION SHIELD INTERFACE ' y 4X y 

3' TMAX =' yFl3.6y 'DEG.F' ) 

217 FORMAT ( ///23X y ' TEMPERATURE PROI [LE IN DEGREES F ') 

215 F0RMAT(///3SXy ' TInE TAU ^='y|IS.6y' SECS') 

216 F0RMAT(/lSXySF13.'5) 

99V STOP 

END 

SUBROUT I NE TRI DAG ( I F y L y A y B y C y D » 0 ) 

DIMENSION A( 1 ) y B< I ) yC( 1 ) ylKl ) yO(l ) yBETAUGO) y GAMMA ( 100) 

DO 10 1=1 y 100 
BETA< I )=0. 

GAMMA ( I )=0. 

10 CONTINUE 

BETA( 1F)=--B( IF) 

GAMMA (IF )=D< 1. F ) / BETA ( IF ) 

IFPl =IF+1 
DO L I = IFPlyl. 
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101 FORMAT <54X» I3/54X» I3/54X»F12. 6) 

102 FORMAK 7X» I2f7X,F10.5»9X»F13.6» 14X»F13.6) 

103 FORMAT (53X»I3/53XfI3/53X>I3/53X. 13) 

104 F0RMAT(53X,I3/53Xf I3/S3X>I3/53X».l3) 

105 F0RMAT(53X»F11 .5/53X»Fll ,5/53X»Fll ,5/53Xf Fll .5) 

106 FORMAT (53X»I3/53Xf 13) 

116 F0RMAT(46X*F13.6/46X»F13.6) 

107 FORMAT ( 46X rF 13. 6/46X*F 15. 6/46X»F13.6/46X»F 15,6) 

100 FORMAT (46XrF13.6/46X>F13. 6) 

109 FORMAT ( 61 X» 1 3/47X r F 1 3 . 6/ 47X » FI 3 . 6/47X r F 1 3 . 6/47X » F 13 . 6/47X » FI 3 . 6 ) 

110 FORMAT ( 4 7X » F 1 3 . 6/ 4 7X f F 1 3 . 6 ) 

111 FORMA r < 50X . F 1 3 . 6/50X » 1 3/50X i- F 1 3 . 6/50X f 1 3 ) 

112 FORMA r ( 50X , F 1 3 . 6/50X . 1 3/ 50X r F 1 3 . 6 ) 

200 FORMAK/ ' fOlAL NUMDER OF SLAHO ' , 31X . ' N = ' k 13/ ' TOTAL NUMBER 
1 OF N0DF3' »31X, 'M=' r 13/ ' TOTAL LENGTH OF OOMT'OGI TE GLAB'y 
219X, ' TLEN=' yF13.6y ' INCUS' ) 

201 FGRMAT(' INTERNAL HEAT GENERATION IN SLAB NUF .'LR ' y 9X » ' IJ='y 

lI3/33Xy ' BETWEEN NOBE N01 = ' t I3/36X y ' AND NODt. N02='» 13) 

202 F0RMAT(//10Xy 'THERE IS NO HEAT SOURCE PRESENT ') 

203 FORMATC/ ' POINT HEAT SOURCE IS PRESENT A T ' y 1 7X y ' NODE^= ' y 13/ 

133Xy ' BETWEEN SLAB N1 = ' y 1 3/36X y ' AND SLAB N2^='yI3) 

204 FORMAT ( / ' CONS ( ANT HE A T INPUT Ol- ' y 28X y ' 0 = ' y f 1 3 . 6 y ' WATTS/ I N# I N ' ) 

224 FORMAT(/ ' ON-TIME FOR HEAT INPU I ' y 27X y ' TON= ' y F13 . 6 y ' SECS ' / 

1' OFF-TIME FUR HEAT INPU T ' y 25X y ' TOFF--= ' y F 1 3 . 6 y ' SECS ' / 

1' VARIABLE HEAT INPUT OF ' y 20X y ' 00- ' y F 1 3 . 6 y ' WA I TS/ IN)K IN ' ) 

205 FORMAK// ' CONSTANT TEMPERAlURE A( X-O'yl/Xy' TX0-'yF13.6y 
1' DEG.F') 

206 FORMAK// ' CONOECTION OCCURS AT X-0 ' /I IX y ' AMBIENT TEMPERAT 

lURE' y5Xy ' TC1-' .F13.6y ' DEG.F' /I IX y ' HEAT TRANSFER COEFF. ' y5Xy 'HI-' y 
2F 15 . 6 y ' B . T . U/TIR . F T . F T . DEG . F ' ) 

207 FORMATC// ' CONSTANT TEMPERATURE AT X- i ' y 1 /X y ' TX1= ' y FI 3 . 6 y 
1 'DEG.F' ) 

208 FORMATC// ' CONOECTIGN OCCURS AT X- K / 1 IX y ' AMBIENT TEMPERAT 

lURE' ySXy 'TG2-' yF13.6y ' DEG.F '/I IX y 'HEAT TRANSFER COEFF. ' ySXy 'H2-' y 
2F15.6y 'B.T.U/HR.FT.FT.DEG.F' ) 

209 FORMATC/ ' THE INITIAL TEMPERATURE IN THE COMPOSITE SLAB TIN -'y 


>? 



BETA<I)=B<I)-A(I)#C<I--1)/BETA<I-1) 

1 GAMMA<I)=<D(I)-ACI)»GAMMA(I-1) )/BETA(I) 
V(L)=GAMHA(L) 

LAST=L“IF 
DO 2 K= If LAST 
I=L-K 

2 V<I)=CAHMA<I)-C(I)^'J<IH)/BETA(I) 

RETURN 

END ' 

FUNCTION QHEATC TAUf TON » TUFF T«V) 
IN=IFIX<TAU/(TQN+TOFF) ) 

IP=IN+1 

A =IN*<TONFTOFF) 

B =IN*TOFFflP»TUN 
C =IP»(TONTTOFF) . ' 

OHEAT=Q 

IF< (B.LT.TAU) ♦ AND. d AU.LT.C) )QHEAT=0. 

IF < TAU . EQ . TON ) QHEAT =QV 

TTON=TON+TOFF 

IF ( TAU . EQ ♦ TTON ) QHEAT=QU 

RETURN 

END 



1. Report No. 


2. Government Accenion No. 


3. Recipient's Cetalog No. 


NASA CR- 165607 


4. Title and Subtitle 

NUMERICAL SIMULATION OF ONE-DIMENSIONAL HEAT 
TRANSFER IN COMPOSITE BODIES WITH PHASE CHANGE 


7. Author(s) 

Kenneth J. DeWitt sind Gurudutt Baliga 


t. Rerforniing Orgenizetion Neme end Addreti 
University of Toledo 
Department of Chemical Engineering 
Toledo, Ohio 43606 


12. Sponsoring Agency Neme end Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 


5. Report Date 

March 1982 


6. Performing Organizetion Code 


I. Performing Organization Report No. 

None 


10. Work Unit No. 


11. Contract or Grant No. 

NAG 3-72 


13. Type of Report end Period Covered 
Contractor Report 


14. Sponsoring Agency Code 

505-44-12 


IS. Suppiementary Notes 

Final report. Project Manager, Robert J. Shaw, Propulsion Systems Division, NASA Lewis Re- 
search Center, Cleveland, Ohio 44135. Based on thesis submitted by Gurudutt Baliga as partial 
fulfillment of the requirements for the degree Master of Science to the University of Toledo, 
Toledo. Ohio in 1980. 


16. Abstract 

A numerical simulation is developed to investigate the one-dimensional heat transfer occurring 
in a system composed of a layered aircraft blade having an ice deposit on its surface. The finite- 
difference representation of the heat conduction equations is done using the Crank-Nicolson im- 
plicit finite-difference formulation. The simulation considers uniform or time-dependent heat 
sources, from heaters which can be either point sources or of finite thickness. For the icewater 
phase change, a numerical method which approximates the latent heat effect by a large heat capac- 
ity over a small temperature interval has been applied. The simulation describes the temperature 
profiles within the various layers of the de-icer pad, as well as the movement of the ice-water 
interface. The simulation could also be used to predict the one-dimensional temperature profiles 
in any composite slab having different boundary conditions. 


17. Kty Words (Suggttttd by Author(s)) 

Aircraft icing 

Electrothermal deicing 

One-dimensional, transient heat conduction 


It. Distribution Stattirwnt 

Unclassified - unlimited 
STAR Category 01 


20. Sacurity Oissif. (of this pige) 

21. No. of Pages 

Unclassified 

88 


It. Stcurity Qissif. (of this rtport) 
Unclassified 


* For sale by the National Technical Information Service, Springfielii, Virginia 22161 

























3 1176 00512 8484 


till) 


t 

i 

* 


I 


i 



